Regresión lineal Datos Vino

Guillermo Villarino

Curso 2022-2023

Preliminares

En este documento ajustaremos algunos modelos de regresión lineal a los datos sobre venta de vinos. Para ello, utilizamos el conjunto de datos que generamos tras la depuración, asegurando un conjunto de datos “limpios” y exentos de ciertos peligros.

En primer lugar leemos el archivo generado en la depuración. Nótese que aquí se podría “repetir” el proceso con distintas posibilidades de tratamiento previo de las variables. Tal vez probar con winsorize e imputaciones aleatorias y luego comparar con otro esquema como el de eliminición de outliers y missing o con conversión de outliers a missings e imputación por modelos multivariantes…Todo un mundo de posibilidades!!

import pandas as pd
import numpy as np
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression

# Leer datos depurados datosvinoDep
vinosDep = pd.read_csv('C:\\Users\\Guille\\Documents\\MineriaDatos_2022_23\\PARTE I_Depuracion y Regresiones\\Dia1_MDDepuracion\\DatosVinoDep_winsRand.csv', index_col=0)

# Descriptivo de comprobación
vinosDep.head()
##    ID  Acidez  AcidoCitrico  Azucar  ...  Region  prop_missings  Beneficio  Compra
## 0   2    0.16         -0.81   26.10  ...     1.0       7.692308        515       1
## 1   4    2.64         -0.88   14.80  ...     3.0       0.000000        585       1
## 2   8    0.29         -0.40   21.50  ...     1.0       0.000000          0       0
## 3  11   -1.22          0.34    1.40  ...     2.0       7.692308        775       1
## 4  12    0.27          1.05   11.25  ...     2.0       0.000000        596       1
## 
## [5 rows x 17 columns]

Solicitamos información sobre el tipo de variables del archivo a modo de comprobación y observamos que se han cambiado los tipos! Esto se produce en el proceso de escritura y lectura del csv… Habrá que cambiarlo.

vinosDep.info()
## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 6365 entries, 0 to 6364
## Data columns (total 17 columns):
##  #   Column          Non-Null Count  Dtype  
## ---  ------          --------------  -----  
##  0   ID              6365 non-null   int64  
##  1   Acidez          6365 non-null   float64
##  2   AcidoCitrico    6365 non-null   float64
##  3   Azucar          6365 non-null   float64
##  4   CloruroSodico   6365 non-null   float64
##  5   Densidad        6365 non-null   float64
##  6   pH              6365 non-null   float64
##  7   Sulfatos        6365 non-null   float64
##  8   Alcohol         6365 non-null   float64
##  9   CalifProductor  6365 non-null   int64  
##  10  PrecioBotella   6365 non-null   float64
##  11  Etiqueta        6365 non-null   object 
##  12  Clasificacion   6365 non-null   object 
##  13  Region          6365 non-null   float64
##  14  prop_missings   6365 non-null   float64
##  15  Beneficio       6365 non-null   int64  
##  16  Compra          6365 non-null   int64  
## dtypes: float64(11), int64(4), object(2)
## memory usage: 895.1+ KB
# Descripción de los datos
vinosDep.describe(include='all')
##                   ID       Acidez  ...    Beneficio       Compra
## count    6365.000000  6365.000000  ...  6365.000000  6365.000000
## unique           NaN          NaN  ...          NaN          NaN
## top              NaN          NaN  ...          NaN          NaN
## freq             NaN          NaN  ...          NaN          NaN
## mean     8010.702278     0.330902  ...   452.380204     0.785232
## std      4654.939139     0.768507  ...   308.380542     0.410694
## min         2.000000    -2.050000  ...     0.000000     0.000000
## 25%      3980.000000     0.130000  ...   236.000000     1.000000
## 50%      8065.000000     0.280000  ...   480.000000     1.000000
## 75%     12027.000000     0.650000  ...   671.000000     1.000000
## max     16128.000000     2.680000  ...  1568.000000     1.000000
## 
## [11 rows x 17 columns]
# Lista de columnas con menos de 10 valores distintos. Potenciales factores!
to_factor = list(vinosDep.loc[:,vinosDep.nunique() < 10]);  

# Podemos cambiar el tipo de todas ellas a factor de una vez
vinosDep[to_factor] = vinosDep[to_factor].astype('category')

# Ordenmaos categorías de los fatcores de interés
vinosDep["Etiqueta"] = vinosDep["Etiqueta"].cat.reorder_categories(['MM','M','R','B','MB'])
vinosDep["Clasificacion"] = vinosDep["Clasificacion"].cat.reorder_categories(['Desc','*','**','***','****'])

Estudio descriptivo de relaciones con la respuesta

En este apartado intentaremos descubrir a priori las relaciones marginales de las variables con la variable objetivo para hacernos una idea de cuales de ellas serán potencialmente influyentes en los modelos de regresión que ajustemos.

import seaborn as sns

# g = sns.PairGrid(vinosDep.iloc[:,6:])
# g.map_diag(sns.histplot)
# g.map_lower(sns.kdeplot)
# g.map_upper(sns.scatterplot)
# g.add_legend()
# 
# del g

Podemos generar este tipo de gráficos de rejilla, donde se pueden visualizar las relaciones a pares de todas las variables del archivo. En este caso, dado el caracter de los gráficos, solo será de aplicación a variables de tipo numérico.

No tiene muy buena pinta ya que las variables no parecen generar un patrón especial frente a la varuiable objetivo…se intuye poco valor predictivo de las variables para la regresión lineal.

Por otro lado, llama la antención la gran cantidad de 0 que tiene la variable objetivo continua Beneficio…Distribución muy compleja de modelizar..Tal vez habría que replantear el análisis y generar un modelo de regresión lineal para estimar el beneficio de los vinos que presentan algún beneficio.. Teniendo así una distribución de la variable respuesta más llevadera.

En este punto, podríamos pensar en una doble vía de actuación para el modelo predictivo final que generemos. Por una parte, un modelo de clasificación previa que estime la probabilidad de ser un vino con Compra = 1 (todo vino que tiene beneficio > 0) y posterioremente, para los clasificados como 1 en compra, un modelo de regresión lineal que estime el beneficio esperado. De alguna forma, estamos planteando un ensamble de modelos. Solo por plantear alternativas.

Variables de control

No es mala idea generar un par de variables de “control” para la evaluación de los efectos de los predictores frente a la respuesta. La idea es la siguiente: si generamos variables en el más estricto sentido aleatorio (por ejemplo siguiendo una distribución uniforme[0,1]) cualquier relación que estas presenten con la variable respuesta serán debidas puramente al azar, con lo que se pueden considerar relaciones espurias, es decir, falsas.

Por tanto, ya sea en la inspección preliminar de relaciones con la respuesta mediante correlación (relación lineal, válido para continua-continua) o VCramer (asociación en tablas de contingencia, válido para cruce de variables categóricas/nominales o continuas tramificadas) o bien en los propios modelos de regresión, las variables que presenten una menor relación con la respuesta que las variables de control, tendrán una sombra de sospecha sobre la veracidad de esa relación y probablemente serán descartadas, al menos en su estado original (siempre se pueden tratar de transformar, tramificar etc)

vinosDep['aleatorio'] = np.random.uniform(0,1,size=vinosDep.shape[0])
vinosDep['aleatorio2'] = np.random.uniform(0,1,size=vinosDep.shape[0])
vinosDep.head()
##    ID  Acidez  AcidoCitrico  Azucar  ...  Beneficio  Compra  aleatorio  aleatorio2
## 0   2    0.16         -0.81   26.10  ...        515       1   0.097808    0.593139
## 1   4    2.64         -0.88   14.80  ...        585       1   0.366505    0.062341
## 2   8    0.29         -0.40   21.50  ...          0       0   0.502511    0.794030
## 3  11   -1.22          0.34    1.40  ...        775       1   0.928091    0.538947
## 4  12    0.27          1.05   11.25  ...        596       1   0.684028    0.511079
## 
## [5 rows x 19 columns]

Interesante función para obetener un reporte descriptivo del archivo completo que evalúa a nivel univariante y bivariante las variables presentes en el archivo. También da información sobre valores extremos y missings, por lo que la podríamos utilizar en nuestra etapa de depuración de los datos.

# pip install pandas_profiling
import pandas_profiling
#pandas_profiling.ProfileReport(vinosDep)

Muchísima información disponible en este reporte automático. A nivel univariante podemos comprobar distribuciones de variables continuas y representatividad de las categorías de los factores. A nivel bivariado, se evalúan las relaciones a pares de distintas formas. Especialmente interesante el gráfico de correlación Auto por su filosofía. Genera la correlación de spearman para continua-continua y v de cramer para categórica-categórica y para continua-categórica previamente tramificada la primera. Va muy en nuestra línea de pensamiento.

Con ese 21% de valores nulos, vamos a decidir partir el archivo para considerar solamente los vinos de Compra = 1 y ajustaremos el modelo de regresión sobre esta submuestra. Es importante que estas decisiones sean consensuadas y se especifiquen claramente en el informe de análisis, trasladándose a la fase de interpretación de resultados.

Creamos el archivo vinosCompra a este efecto.

vinosCompra = vinosDep.loc[vinosDep.Compra == 1].drop(['Compra'],axis=1)
vinosCompra.info()
## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 18 columns):
##  #   Column          Non-Null Count  Dtype   
## ---  ------          --------------  -----   
##  0   ID              4998 non-null   int64   
##  1   Acidez          4998 non-null   float64 
##  2   AcidoCitrico    4998 non-null   float64 
##  3   Azucar          4998 non-null   float64 
##  4   CloruroSodico   4998 non-null   float64 
##  5   Densidad        4998 non-null   float64 
##  6   pH              4998 non-null   float64 
##  7   Sulfatos        4998 non-null   float64 
##  8   Alcohol         4998 non-null   float64 
##  9   CalifProductor  4998 non-null   int64   
##  10  PrecioBotella   4998 non-null   float64 
##  11  Etiqueta        4998 non-null   category
##  12  Clasificacion   4998 non-null   category
##  13  Region          4998 non-null   category
##  14  prop_missings   4998 non-null   category
##  15  Beneficio       4998 non-null   int64   
##  16  aleatorio       4998 non-null   float64 
##  17  aleatorio2      4998 non-null   float64 
## dtypes: category(4), float64(11), int64(3)
## memory usage: 606.0 KB
#pandas_profiling.ProfileReport(vinosCompra)

En el reporte para el nuevo archivo, volvemos a tener alta correlación (en un sentido amplio) de las variables Etiqueta y Clasificación con la variable objetivo, apareciendo algo coloreada la variable Alcohol. Pocas relaciones destacables a parte de esto.

Ya que vamos a hacer cosas como evaluación de las relaciones entre los predictores y la respuesta o creación masiva de transformaciones para conseguir linealidad, lo mejor es separar las respuestas y quedarnos con el input depurado, de esta forma podemos aplicar una misma función a todo el conjunto sin peligro de transformar las respuestas y cosas raras que puedan suceder.

# Eliminar variable objetivo continua 
varObjCont = vinosCompra.Beneficio
imputCompra = vinosCompra.drop(['ID','Beneficio'],axis=1)
varObjCont.info()
## <class 'pandas.core.series.Series'>
## Int64Index: 4998 entries, 0 to 6364
## Series name: Beneficio
## Non-Null Count  Dtype
## --------------  -----
## 4998 non-null   int64
## dtypes: int64(1)
## memory usage: 78.1 KB

Nos encantaría tener un procedimiento que evaluara las relaciones de todas las variables del input frente a la variable objetivo (ya sea continua o binaria) y que de alguna forma ordenara esos efectos por intensiadad de asociación para generar un ranking tentativo a priori de las vaiables más interesantes a nivel individual (es importante destacar que estas relaciones solamente cuentan el efecto marginal de una variable con la objetivo sin considerar las posibles interacciones existentes entre variables que pueden modificar en modelo estas asociaciones).

Vamos a definir algunas funciones que nos faciliten el proceso de búsqueda, automatizando lo más posible para que luego sea aplicar y evaluar!

Primero nos generamos la función cramers_v, cuyo objetivo es calcular el valor de v de cramer para la asociación entre dos variables cualesquiera (se entiende que una de ellas será una variable objetivo…pero podría funcionar entre predictores de igual forma..) Molaría que la función distinguiera el caso de variable continua y categórica puesto que la asociación mediante vCramer se entiende en términos de cruce de dos factores y, por tanto, si la variable es continua, debe existir un proceso previo de tramificación o discretización de la variable. Entonces, programamos el proceso de tal manera que si encuentra columnas numéricas las tramifique (en 5 tramos por ejemplo) y sino simplemente calcule el valor cramer que se extrae de la tabla de contingencia de los dos factores introducidos.

La filosofía de programación es la misma que en otras funciones definidas anteriormente. Como realmente querremos aplicarla sobre el input completo, estamos pensando en utilizar un apply para que aplique un mismo proceso a todas las columnas que vaya encontrando, manteniendo la variable objetivo. Por ello, es importante separar esta objetivo del conjunto de predictores así como lo hacíamos en imputaciones y transformación de missings.

# Librería estadística!
import scipy.stats as stats

# Función para calcular VCramer (dos nominales de entrada!)
def cramers_v(var1, varObj):
    
    if not var1.dtypes == 'category':
        bins = min(5,var1.value_counts().count())
        var1 = pd.cut(var1, bins = 5)
    if not varObj.dtypes == 'category': #np.issubdtype(varObj, np.number):
        bins = min(5,varObj.value_counts().count())
        varObj = pd.cut(varObj, bins = 5)
        
    data = pd.crosstab(var1, varObj).values
    vCramer = stats.contingency.association(data, method = 'cramer')
    return vCramer


# Ejemplo uso univariante
#cramers_v(vinosCompra['Etiqueta'],vinosCompra['Beneficio'])

# Aplicar la función al input completo contra la objetivo
tablaCramer = pd.DataFrame(imputCompra.apply(lambda x: cramers_v(x,varObjCont)),columns=['VCramer'])

# Obtener el gráfico de importancia de las variables frente a la objetivo continua según vcramer
import plotly.express as px
px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones frente al Beneficio').update_yaxes(categoryorder="total ascending")

Ya tenemos nuestro ranking de variables según cramer para todos los predictores en relación a la variable objetivo continua, Beneficio. Como ya intuíamos, las variables categóricas Etiqueta y Clasificación son las que mayor relación presentan, es decir, las que mayor patrón de asociación pueden aportar para el futuro modelo.

Entendiendo el test Chi-cuadrado

Vamos a realizar un análisis de asociación de una tabla de contingencia basada en Chi-cuadrado. Para ilustrar esto, nos permitimos el lujo de adelantar algo que observaremos posteriormente y que podría afectar a nuestro modelo por la presencia de interacción entre las variables más relevantes Clasificación y Etiqueta.

# Opción para calcular las medias de Beneficio en los cruces de Clasificación y etiqueta
pd.pivot_table(vinosCompra,index =['Etiqueta'],columns=['Clasificacion'],values=['Beneficio'],aggfunc=[np.mean])

# Tabla de contingencia con la cuenta de frecuencia de observaciones en cada cruce
##                      mean                                                
##                 Beneficio                                                
## Clasificacion        Desc           *          **         ***        ****
## Etiqueta                                                                 
## MM             227.714286  292.464646  349.676471  352.909091         NaN
## M              365.428571  403.221445  449.474178  523.120968  635.846154
## R              507.620553  536.413989  581.393564  636.366733  725.234694
## B              706.515625  635.623457  701.191919  775.459016  857.936416
## MB             900.750000  830.000000  826.568182  857.851351  929.973684
pd.crosstab(vinosCompra.Etiqueta,vinosCompra.Clasificacion)
## Clasificacion  Desc    *   **  ***  ****
## Etiqueta                                
## MM               63   99   34   11     0
## M               273  429  426  124    13
## R               253  529  808  499    98
## B                64  162  396  366   173
## MB               12   10   44   74    38

Vayamos un poco más allá por un momento y nos planteamos si habrá un efecto de interacción entre estas dos variable relevantes. Que demonios queremos decir con esto? Si hay interacción, es posible que el efecto conjunto de las dos variable sobre la respuesta no sea uniforme o bien no sea el resultado de la adición simple de ambos efectos. De esta forma, nos planteamos si es igualmente probable encontrar un vino con Clasificación máxima y la peor Etiqueta, que un vino con Clasificación máxima y la mejor de las Etiquetas… La lógica puede inducir a pensar que debería existir cierto patrón de asociación..

# Probamos el test chi cuadrado para ver las salidas
stats.chi2_contingency(pd.crosstab(vinosCompra.Etiqueta,vinosCompra.Clasificacion).values)

#pd.crosstab(data_trans.Beneficio_cut,vinosCompra.Clasificacion,margins=True,normalize='columns')
## (842.5114902984811, 5.340776470642584e-169, 16, array([[ 27.54201681,  50.90096038,  70.7394958 ,  44.48139256,
##          13.33613445],
##        [168.31232493, 311.06142457, 432.29691877, 271.83073229,
##          81.49859944],
##        [290.98739496, 537.77971188, 747.37815126, 469.95558223,
##         140.89915966],
##        [154.47478992, 285.4879952 , 396.75630252, 249.48259304,
##          74.79831933],
##        [ 23.68347339,  43.76990796,  60.82913165,  38.24969988,
##          11.46778711]]))
  • Primer valor: Valor del estadístico Chi-cuadrado
  • Segundo valor: p-valor del contraste de hipótesis (recordemos H0: no existe asociación)
  • Tercer valor: Grados de libertad de la distribución de contraste (filas-1)*(columnas-1)
  • Cuarto valor: Matriz de valores esperados bajo H0 (según las distribuciones marginales de las variables, cuan probable es encontrar casilla 1,1… p(x=1)*p(y=1)

El valor de chi es muy alto, por lo que el p-valor es muy pequeño, prácticamente 0. Por lo que se evidencia una asociación significativa entre estas dos variables, o bien, la matriz de valores observados es significativamente distinta a la matriz de valores esperados bajo independencia.

Precaución!! Hay una cosa importante en el test Chi-cuadrado. Los resultados no son robustos si alguna de las casillas tiene frecuencia esperada menor que 5. Si se da este caso, sería necesario recodificar alguna de las variables del cruce para obtener tamaño muestral suficiente en las marginales. En nuestro caso, esto no llega a suceder (ver matriz de frecuencias esperadas), sin embargo tenemos una casilla con frecuencia observada 0, lo cual puede estar tirando mucho del valor de la chi y de alguna forma distorsionando sus resultados. Por ello, no es mala práctica unir alguna categoría para evitar esto.

Veamos los cruces de estas categóricas con Beneficio.

#sns.boxplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
print(sns.violinplot(x='Etiqueta',y='Beneficio',data=vinosCompra,palette='viridis'))
#sns.stripplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
#sns.swarmplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
## AxesSubplot(0.125,0.11;0.775x0.77)
#sns.boxplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')
#sns.violinplot(x='Etiqueta',y='Beneficio',data=vinosCompra,palette='viridis')
sns.stripplot(x='Clasificacion',y='Beneficio',data=vinosDep,palette='viridis')
#sns.swarmplot(x='Etiqueta',y='Beneficio',data=vinosDep,palette='viridis')

Evaluación de la asociación categorica-numérica mediante ANOVA

Los contrastes estadísticos pueden resultar de utilidad para contatar la significación de las relaciones visualizadas a nivel gráfico. De esta forma, si tenemos una factor y queremos cuantificar esta significación con respecto a la variable objetivo continua, no es descabellado pensar en un análisis de la varianza con un factor para valorar si las distribuciones en los subgrupos difieren de la distribución global.

El Anova realiza la descomposición de la suma de cuadrados en las componentes explicada y residual, de tal forma que imputamos una parte de la variabilidad de la variable continua al cambio en los niveles del factor y la parte que no se puede explicar solamente por ese cambio, queda como variabilidad no explicada. La idea es que la variabilidad explicada sea mucho mayor que la residual.

Para la variable Región: Se observa que el p-valor es superior a 0,05 por lo que la asociación no es significativa al 95% de nivel de confianza. Esta variable no influye mucho en el Beneficio.

# get ANOVA table as R like output
import statsmodels.api as sm
from statsmodels.formula.api import ols

# ANOVA región
model = ols('Beneficio ~ C(Region)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table
##                  sum_sq      df         F    PR(>F)
## C(Region)  2.362961e+05     2.0  2.372917  0.093313
## Residual   2.487021e+08  4995.0       NaN       NaN

Para la variable Etiqueta: Se observa que el p-valor es 0 por lo que la asociación es significativa a cualquier nivel de confianza. Esta variable influye mucho en el Beneficio.

# ANOVA Etiqueta
model = ols('Beneficio ~ C(Etiqueta)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table
##                    sum_sq      df           F  PR(>F)
## C(Etiqueta)  9.240747e+07     4.0  736.899964     0.0
## Residual     1.565309e+08  4993.0         NaN     NaN

ANOVA de doble vía. Puede resultar útil para la evaluación de efectos conjuntos de factores sobre la objetivo continua. En este caso, evaluamos la interacción de la que ya sospechábamos entre Clasificación y Etiqueta.

# ANOVA con dos factores e interacción
model = ols('Beneficio ~ C(Etiqueta) + C(Clasificacion) + C(Etiqueta):C(Clasificacion)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table
##                                     sum_sq      df           F        PR(>F)
## C(Etiqueta)                   7.500830e+07     4.0  659.687718  0.000000e+00
## C(Clasificacion)              4.971909e+06     4.0   43.727257  4.012464e-36
## C(Etiqueta):C(Clasificacion)  1.016574e+07    16.0   22.351569  7.635523e-64
## Residual                      1.413894e+08  4974.0         NaN           NaN

Se constata la significación estadística de todos los efectos presentes en el modelo. Podríamos decir que existe influencia marginal de ambos factores sobre el Beneficio pero también una influencia conjunta de ambos! Así, el beneficio esperado puede depender de la combinación de Clasificación y Etiqueta del vino.

Para ilustrar este efecto, podemos recurrir al gráfico de interacciones entre dos factores y la variable continua. Aquí, detalle de programación y es que no acepta el tipo category (con su orden) por lo que habrá que convertir las columnas en tipo ‘object’ (pero solo para esto!)

from statsmodels.graphics.factorplots import interaction_plot
import matplotlib.pyplot as plt


# Gráfico de interacción de las componentes del ANOVA
fig = interaction_plot(x=vinosCompra.Etiqueta.astype('object'), trace=vinosCompra['Clasificacion'].astype('object'), response=vinosCompra['Beneficio'])
plt.show()

Se evidencia interacción entre los factores cuando las líneas se cortan en algún punto. Esto quier decir que la evolución del Beneficio a lo largo de las Etiquetas no se produce de la misma forma en todas los niveles de Clasificación.

Se puede observar esa falta de dato en la combinación MM - ****

Vamos a hacer la prueba de recategorizar alguno de los niveles implicados. Para decidir cual de ellos, recurriremos a 1) frecuencia de la categoría (reporte) y 2) Similitud con la categoría adyacente (boxplot-violin). De esta forma decidimos que parece más conveniente unir Etiqueta en sus categorías M y MM.

# Unimos la categoría minoritaria MM con M en Etiqueta
vinosCompra['Etiqueta_r'] = vinosCompra.Etiqueta.replace(['M','MM'],'M-MM',inplace=False)
pd.crosstab(vinosCompra.Etiqueta_r,vinosCompra.Clasificacion)
## Clasificacion  Desc    *   **  ***  ****
## Etiqueta_r                              
## M-MM            336  528  460  135    13
## R               253  529  808  499    98
## B                64  162  396  366   173
## MB               12   10   44   74    38
# Probamos el test chi cuadrado para ver las salidas
stats.chi2_contingency(pd.crosstab(vinosCompra.Etiqueta_r,vinosCompra.Clasificacion).values)
## (800.5770302943753, 1.244618385605778e-163, 12, array([[195.85434174, 361.96238495, 503.03641457, 316.31212485,
##          94.83473389],
##        [290.98739496, 537.77971188, 747.37815126, 469.95558223,
##         140.89915966],
##        [154.47478992, 285.4879952 , 396.75630252, 249.48259304,
##          74.79831933],
##        [ 23.68347339,  43.76990796,  60.82913165,  38.24969988,
##          11.46778711]]))
# Gráfico de interacción de las componentes del ANOVA
fig = interaction_plot(x=vinosCompra.Etiqueta_r.astype('object'), trace=vinosCompra['Clasificacion'].astype('object'), response=vinosCompra['Beneficio'])
plt.show()

Sigue evidenciándose interacción entre los factores, especialmente en la categoría desconocido.

Transformaciones de las variables continuas

El principal objetivo de las transformaciones de las variables continuas es conseguir linealidad frente a la variable objetivo. De esta forma se prueban las transformaciones típicas (log, exp, potencias y raíces) y se escoge aquella que mayor coeficiente de correlación o valor V de Cramer presenta con la respuesta continua (o binaria en su caso).

Definamos pues una función que nos facilite este proceso, de tal manera que aplique sobre una columna genérica y una variable objetivo e implemente la posibilidad de actuar en base al coeficiente de correlación (relación lineal de continua-continua) o al valor V de cramer (asociación general entre categórica-categórcia).

La función va a hacer una traslación a valores positivos de la variable que encuentre (predictor continuo siempre!! no tiene sentido el log de una categórica!) y seguidamente va a generar las posibles transfromaciones típicas (log, exp, potencias y raíces 2 y 4) para aplicar estos criterios sobre todas las posibilidades y escoger la que mejor valor obtenga. Vamos a ello!

Aquí hay un walkaround para conseguir que los nombres de las transfromadas se conserven en el dataset de salida que es fruto de un apply al input continuo. Seguro se puede mejorar!

# Busca la transformación de variables input de intervalo que maximiza la VCramer o la correlación tipo Pearson con la objetivo

def mejorTransf (vv,target, name=False, tipo = 'cramer', graf=False):
    # Traslación a valores positivos de la variable (sino falla log y las raíces!)
    vv = vv + abs(min(vv))+0.0001
      
    # Definimos y calculamos las tranformacione típicas  
    transf = pd.DataFrame({vv.name + '_ident': vv, vv.name + '_log': np.log(vv), vv.name + '_exp': np.exp(vv), 
                         vv.name + '_sqr': np.square(vv), vv.name + '_cuarta': vv**4, vv.name + '_raiz4': vv**(1/4)})
      
    # Distinguimos caso cramer o caso correlación
    if tipo == 'cramer':
      # Aplicar la función cramers_v a cada trasnformación frente a la respuesta
      tablaCramer = pd.DataFrame(transf.apply(lambda x: cramers_v(x,target)),columns=['VCramer'])
      
      # Si queremos gráfico, muestra comparativa entre las posibilidades
      if graf: px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones frente a ' + target.name).update_yaxes(categoryorder="total ascending").show()
      # Identificar mejor transfromación
      best = tablaCramer.query('VCramer == VCramer.max()').index
      ser = transf[best].squeeze()
    
    if tipo == 'cor':
      # Aplicar coeficiente de correlacion a cada trasnformación frente a la respuesta
      tablaCorr = pd.DataFrame(transf.apply(lambda x: np.corrcoef(x,target)[0,1]),columns=['Corr'])
      # Si queremos gráfico, muestra comparativa entre las posibilidades
      if graf: px.bar(tablaCorr,x=tablaCorr.Corr,title='Relaciones frente a ' + target.name).update_yaxes(categoryorder="total ascending").show()
      # identificar mejor transfromación
      best = tablaCorr.query('Corr.abs() == Corr.abs().max()').index
      ser = transf[best].squeeze()
  
    # Aquí distingue si se devuelve la variable transfromada o solamente el nombre de la transfromacion
    if name:
      return(ser.name)
    else:
      return(ser)

# Ejemplo de uso univariante
tr = mejorTransf(vinosCompra.Azucar,varObjCont, tipo='cor')

tr.head(),vinosCompra.Azucar.head()
## (0    15055.31454
## 1    12409.98228
## 3     9604.01960
## 4    11631.64407
## 5     9682.57968
## Name: Azucar_sqr, dtype: float64, 0    26.10
## 1    14.80
## 3     1.40
## 4    11.25
## 5     1.80
## Name: Azucar, dtype: float64)
#tr=pd.DataFrame()
#for i in range(len(imputCompra.select_dtypes(include=np.number))):
#    tr.iloc[:,i] = mejorTransfVCramer(imputCompra[column],varObjCont)
imputCompra.info()
## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 16 columns):
##  #   Column          Non-Null Count  Dtype   
## ---  ------          --------------  -----   
##  0   Acidez          4998 non-null   float64 
##  1   AcidoCitrico    4998 non-null   float64 
##  2   Azucar          4998 non-null   float64 
##  3   CloruroSodico   4998 non-null   float64 
##  4   Densidad        4998 non-null   float64 
##  5   pH              4998 non-null   float64 
##  6   Sulfatos        4998 non-null   float64 
##  7   Alcohol         4998 non-null   float64 
##  8   CalifProductor  4998 non-null   int64   
##  9   PrecioBotella   4998 non-null   float64 
##  10  Etiqueta        4998 non-null   category
##  11  Clasificacion   4998 non-null   category
##  12  Region          4998 non-null   category
##  13  prop_missings   4998 non-null   category
##  14  aleatorio       4998 non-null   float64 
##  15  aleatorio2      4998 non-null   float64 
## dtypes: category(4), float64(11), int64(1)
## memory usage: 656.9 KB
# Aplicar a las variables continuas la mejor transfromación según correlacion frente a varObjCont
transf_cor = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont, tipo='cor'))
# Pedir los nombres de las transromadas
transf_cor_names = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont,tipo='cor', name=True))
# Asignar nombres a las columnas de salida del proceso
transf_cor.columns = transf_cor_names.values
transf_cor
##       Acidez_sqr  AcidoCitrico_exp  ...  aleatorio_log  aleatorio2_cuarta
## 0       4.884542          4.179117  ...      -2.323380           0.123920
## 1      21.997038          3.896583  ...      -1.003377           0.000015
## 3       0.689066         13.198458  ...      -0.074482           0.084480
## 4       5.382864         26.845548  ...      -0.379560           0.068321
## 5       3.349266         13.875157  ...      -2.186242           0.007613
## ...          ...               ...  ...            ...                ...
## 6357    4.973346         13.198458  ...      -0.119273           0.018796
## 6361    7.508148         10.278969  ...      -2.864621           0.004033
## 6362    5.664876          3.669664  ...      -1.735340           0.002771
## 6363    4.928844          3.127081  ...      -0.338840           0.207337
## 6364    5.617374         10.592011  ...      -3.952220           0.086490
## 
## [4998 rows x 12 columns]
transf_cor
##       Acidez_sqr  AcidoCitrico_exp  ...  aleatorio_log  aleatorio2_cuarta
## 0       4.884542          4.179117  ...      -2.323380           0.123920
## 1      21.997038          3.896583  ...      -1.003377           0.000015
## 3       0.689066         13.198458  ...      -0.074482           0.084480
## 4       5.382864         26.845548  ...      -0.379560           0.068321
## 5       3.349266         13.875157  ...      -2.186242           0.007613
## ...          ...               ...  ...            ...                ...
## 6357    4.973346         13.198458  ...      -0.119273           0.018796
## 6361    7.508148         10.278969  ...      -2.864621           0.004033
## 6362    5.664876          3.669664  ...      -1.735340           0.002771
## 6363    4.928844          3.127081  ...      -0.338840           0.207337
## 6364    5.617374         10.592011  ...      -3.952220           0.086490
## 
## [4998 rows x 12 columns]
# Aplicar a las variables continuas la mejor transfromación según cramer frente a varObjCont
transf_cramer = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont, tipo='cramer'))
transf_cramer_names = imputCompra.select_dtypes(include=np.number).apply(lambda x: mejorTransf(x,varObjCont,tipo='cramer', name=True))
transf_cramer.columns = transf_cramer_names.values
transf_cramer
##       Acidez_sqr  AcidoCitrico_cuarta  ...  aleatorio_cuarta  aleatorio2_cuarta
## 0       4.884542             4.182786  ...      9.201856e-05           0.123920
## 1      21.997038             3.422026  ...      1.806986e-02           0.000015
## 3       0.689066            44.314531  ...      7.423561e-01           0.084480
## 4       5.382864           117.175386  ...      2.190970e-01           0.068321
## 5       3.349266            47.850783  ...      1.592608e-04           0.007613
## ...          ...                  ...  ...               ...                ...
## 6357    4.973346            44.314531  ...      6.205854e-01           0.018796
## 6361    7.508148            29.478015  ...      1.055951e-05           0.004033
## 6362    5.664876             2.856979  ...      9.669522e-04           0.002771
## 6363    4.928844             1.689553  ...      2.578546e-01           0.207337
## 6364    5.617374            31.025702  ...      1.362358e-07           0.086490
## 
## [4998 rows x 12 columns]
# Input con transfromadas según correlación
imput_transf_cor = imputCompra.join(transf_cor)

# Input con transfromadas según cramer
imput_transf_cramer = imputCompra.join(transf_cramer)
imput_transf_cramer.info()
## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 4998 entries, 0 to 6364
## Data columns (total 28 columns):
##  #   Column                Non-Null Count  Dtype   
## ---  ------                --------------  -----   
##  0   Acidez                4998 non-null   float64 
##  1   AcidoCitrico          4998 non-null   float64 
##  2   Azucar                4998 non-null   float64 
##  3   CloruroSodico         4998 non-null   float64 
##  4   Densidad              4998 non-null   float64 
##  5   pH                    4998 non-null   float64 
##  6   Sulfatos              4998 non-null   float64 
##  7   Alcohol               4998 non-null   float64 
##  8   CalifProductor        4998 non-null   int64   
##  9   PrecioBotella         4998 non-null   float64 
##  10  Etiqueta              4998 non-null   category
##  11  Clasificacion         4998 non-null   category
##  12  Region                4998 non-null   category
##  13  prop_missings         4998 non-null   category
##  14  aleatorio             4998 non-null   float64 
##  15  aleatorio2            4998 non-null   float64 
##  16  Acidez_sqr            4998 non-null   float64 
##  17  AcidoCitrico_cuarta   4998 non-null   float64 
##  18  Azucar_raiz4          4998 non-null   float64 
##  19  CloruroSodico_cuarta  4998 non-null   float64 
##  20  Densidad_log          4998 non-null   float64 
##  21  pH_cuarta             4998 non-null   float64 
##  22  Sulfatos_sqr          4998 non-null   float64 
##  23  Alcohol_ident         4998 non-null   float64 
##  24  CalifProductor_log    4998 non-null   float64 
##  25  PrecioBotella_sqr     4998 non-null   float64 
##  26  aleatorio_cuarta      4998 non-null   float64 
##  27  aleatorio2_cuarta     4998 non-null   float64 
## dtypes: category(4), float64(23), int64(1)
## memory usage: 1.1 MB

Podemos comprobar si las transformaciones han aumentado el valor de correlación lineal con la respuesta.

En el gráfico de correlaciones se intuye la baja relación lineal presente entre las variables del archivo 2 a 2. Esto nos hace pensar que las variables continuas, a diferencia de las categóricas, no van a tener demasiada influencia en los modelos de regresión frente a la respuesta (véase el poco color que presenta la fila 1). Por otro lado, podemos estar respirar con cierta tranquilidad ante el hecho de la baja relación entre los predictores (véase la ausencia de color en todo el gráfico), cosa que evitará problemas de colinealidad en los modelos.

Es evidente que, las relaciones de cada variable con su transformada han de ser muy elevadas, ya que hay un mecanismo que genera una en función exacta de la otra, por lo que esto es normal.

Nota: La principal precaución que hay que tener es no considerar un modelo completo con todas al mismo tiempo puesto que se pueden generar los problemas de colinealidad. Solamente utilizaremos el set completo de variables cuando hagamos un proceso de selección automática de variables, proceso en el cual se elegirán las que más R2 aporten al modelo.

# Matriz de correlaciones incluyendo transformaciones
corr = imput_transf_cor.join(varObjCont).corr()
corr.style.background_gradient(cmap='coolwarm').format(precision=3)
#sns.heatmap(corr, annot=True)
  Acidez AcidoCitrico Azucar CloruroSodico Densidad pH Sulfatos Alcohol CalifProductor PrecioBotella aleatorio aleatorio2 Acidez_sqr AcidoCitrico_exp Azucar_sqr CloruroSodico_log Densidad_log pH_log Sulfatos_exp Alcohol_raiz4 CalifProductor_raiz4 PrecioBotella_sqr aleatorio_log aleatorio2_cuarta Beneficio
Acidez 1.000 -0.023 -0.022 0.006 0.007 0.015 0.027 -0.003 0.009 -0.002 0.005 -0.008 0.961 -0.014 -0.030 -0.007 0.007 0.017 0.010 -0.007 0.003 -0.004 0.013 -0.004 -0.044
AcidoCitrico -0.023 1.000 0.007 -0.017 -0.014 -0.010 0.002 0.003 0.061 -0.013 0.015 -0.005 -0.029 0.796 0.005 -0.006 -0.014 -0.011 0.011 0.006 0.058 -0.015 0.021 -0.010 0.002
Azucar -0.022 0.007 1.000 -0.010 0.013 0.003 -0.001 -0.028 -0.003 -0.020 -0.011 0.014 -0.023 -0.001 0.956 0.004 0.013 0.002 -0.009 -0.032 0.005 -0.018 -0.016 0.004 -0.013
CloruroSodico 0.006 -0.017 -0.010 1.000 0.030 -0.016 -0.004 -0.019 0.014 -0.007 0.004 0.004 0.005 -0.016 -0.014 0.628 0.030 -0.018 -0.009 -0.018 0.023 -0.007 0.009 -0.001 -0.003
Densidad 0.007 -0.014 0.013 0.030 1.000 -0.016 -0.012 0.000 0.035 0.008 0.008 0.015 0.006 -0.013 0.011 -0.002 1.000 -0.015 -0.012 0.007 0.031 0.010 0.010 0.018 -0.015
pH 0.015 -0.010 0.003 -0.016 -0.016 1.000 0.004 -0.017 -0.078 -0.010 -0.003 0.013 0.015 0.012 0.002 -0.001 -0.016 0.990 0.000 -0.015 -0.080 -0.009 0.000 0.025 0.018
Sulfatos 0.027 0.002 -0.001 -0.004 -0.012 0.004 1.000 -0.003 0.019 0.004 0.006 -0.011 0.023 -0.003 -0.006 0.001 -0.012 0.006 0.698 -0.004 0.015 0.005 0.012 -0.013 -0.005
Alcohol -0.003 0.003 -0.028 -0.019 0.000 -0.017 -0.003 1.000 -0.025 0.015 0.002 0.018 0.001 0.008 -0.028 -0.011 0.001 -0.021 -0.010 0.955 -0.033 0.012 0.005 0.018 0.082
CalifProductor 0.009 0.061 -0.003 0.014 0.035 -0.078 0.019 -0.025 1.000 0.019 0.023 -0.016 0.005 0.048 -0.005 0.002 0.035 -0.076 0.007 -0.018 0.890 0.022 0.024 -0.026 -0.052
PrecioBotella -0.002 -0.013 -0.020 -0.007 0.008 -0.010 0.004 0.015 0.019 1.000 0.007 -0.012 -0.006 -0.018 -0.016 0.015 0.008 -0.012 -0.002 0.010 0.022 0.978 0.015 -0.015 0.011
aleatorio 0.005 0.015 -0.011 0.004 0.008 -0.003 0.006 0.002 0.023 0.007 1.000 0.016 0.004 0.017 -0.008 -0.013 0.008 -0.003 -0.007 0.009 0.013 0.002 0.872 0.018 -0.001
aleatorio2 -0.008 -0.005 0.014 0.004 0.015 0.013 -0.011 0.018 -0.016 -0.012 0.016 1.000 -0.007 0.002 0.014 -0.019 0.016 0.013 -0.015 0.012 -0.015 -0.014 0.021 0.866 0.017
Acidez_sqr 0.961 -0.029 -0.023 0.005 0.006 0.015 0.023 0.001 0.005 -0.006 0.004 -0.007 1.000 -0.018 -0.030 -0.012 0.006 0.017 0.006 -0.004 -0.001 -0.008 0.011 -0.001 -0.044
AcidoCitrico_exp -0.014 0.796 -0.001 -0.016 -0.013 0.012 -0.003 0.008 0.048 -0.018 0.017 0.002 -0.018 1.000 -0.004 -0.005 -0.012 0.012 -0.000 0.011 0.049 -0.017 0.025 -0.005 0.014
Azucar_sqr -0.030 0.005 0.956 -0.014 0.011 0.002 -0.006 -0.028 -0.005 -0.016 -0.008 0.014 -0.030 -0.004 1.000 0.005 0.011 0.001 -0.013 -0.033 0.005 -0.012 -0.012 0.002 -0.013
CloruroSodico_log -0.007 -0.006 0.004 0.628 -0.002 -0.001 0.001 -0.011 0.002 0.015 -0.013 -0.019 -0.012 -0.005 0.005 1.000 -0.002 -0.003 0.004 -0.014 0.014 0.016 -0.003 -0.017 0.009
Densidad_log 0.007 -0.014 0.013 0.030 1.000 -0.016 -0.012 0.001 0.035 0.008 0.008 0.016 0.006 -0.012 0.011 -0.002 1.000 -0.015 -0.012 0.008 0.031 0.010 0.010 0.019 -0.015
pH_log 0.017 -0.011 0.002 -0.018 -0.015 0.990 0.006 -0.021 -0.076 -0.012 -0.003 0.013 0.017 0.012 0.001 -0.003 -0.015 1.000 0.001 -0.017 -0.079 -0.011 -0.001 0.026 0.019
Sulfatos_exp 0.010 0.011 -0.009 -0.009 -0.012 0.000 0.698 -0.010 0.007 -0.002 -0.007 -0.015 0.006 -0.000 -0.013 0.004 -0.012 0.001 1.000 -0.010 0.007 0.000 -0.002 -0.008 -0.009
Alcohol_raiz4 -0.007 0.006 -0.032 -0.018 0.007 -0.015 -0.004 0.955 -0.018 0.010 0.009 0.012 -0.004 0.011 -0.033 -0.014 0.008 -0.017 -0.010 1.000 -0.024 0.007 0.015 0.015 0.083
CalifProductor_raiz4 0.003 0.058 0.005 0.023 0.031 -0.080 0.015 -0.033 0.890 0.022 0.013 -0.015 -0.001 0.049 0.005 0.014 0.031 -0.079 0.007 -0.024 1.000 0.023 0.015 -0.021 -0.053
PrecioBotella_sqr -0.004 -0.015 -0.018 -0.007 0.010 -0.009 0.005 0.012 0.022 0.978 0.002 -0.014 -0.008 -0.017 -0.012 0.016 0.010 -0.011 0.000 0.007 0.023 1.000 0.010 -0.019 0.012
aleatorio_log 0.013 0.021 -0.016 0.009 0.010 0.000 0.012 0.005 0.024 0.015 0.872 0.021 0.011 0.025 -0.012 -0.003 0.010 -0.001 -0.002 0.015 0.015 0.010 1.000 0.027 -0.004
aleatorio2_cuarta -0.004 -0.010 0.004 -0.001 0.018 0.025 -0.013 0.018 -0.026 -0.015 0.018 0.866 -0.001 -0.005 0.002 -0.017 0.019 0.026 -0.008 0.015 -0.021 -0.019 0.027 1.000 0.027
Beneficio -0.044 0.002 -0.013 -0.003 -0.015 0.018 -0.005 0.082 -0.052 0.011 -0.001 0.017 -0.044 0.014 -0.013 0.009 -0.015 0.019 -0.009 0.083 -0.053 0.012 -0.004 0.027 1.000

No parece que haya mejorado mucho…Este estudio nos lleva a pensar que las variables continuas no van a aportar mucho a nuestro modelo de regresión, si bien puede suceder que en presencia de otras variables (por ejemplo las categóricas) la influencia de éstas pueda aumentar debido al distinto comportamiento (pendiente) que presenten con la respuesta en los distintos grupos que forma la variable categórica.

Por este motivo, es positivo complementar este análisis descriptivo con la información que dan los modelos sobre la importancia de variables, ya que dentro de modelo, la influencia de la variable se entiende a niveles constantes de todas las demás presentes en el mismo.

En cualquier caso, una vía que se puede explorar aquí para conseguir mayor influencia de las continuas es intentar tramificarlas (pasarlas a categóricas) de tal manera que se puedan descubrir patrones muy alejados d la linealidad que esta puedan contener. En el mundo de la tramificación existen dos estrategias fundamentales,

  1. tramos ad-hoc para una mejor interpretabilidad (por ejemplo tramos de edad acordes a los que se utilizan en las estadísticas oficiales) en los que normalmente se realizan particiones más o menos equidistantes o basadas en un fundamento anterior.

  2. tramos que maximicen cierta relación con la variable objetivo. Se pueden generar los puntos de corte de la variable a tramificar tales que los grupos creados sean lo más distintos entre sí en distribución de la variable respuesta, con lo que se “aseguran” tramos con capacidad de discriminar frente al objetivo. Suelen funcionar bien este tipo de tramificaciones pero tienen cierto peligro de sobreajustar a los datos de entrenamiento…Una de las formas de hacerlo es mediante árboles de regresión/clasificación.

Un árbol realiza particiones binarias de las variables con la premisa de encontrar las mayores diferencias entre los grupos que va formando, justo lo que buscamos!

# http://gnpalencia.org/optbinning/binning_continuous.html
from optbinning import ContinuousOptimalBinning

optb = ContinuousOptimalBinning(name='Azucar', dtype="numerical", max_n_bins=3)
optb.fit(imput_transf_cor['Azucar'].values, varObjCont)
ContinuousOptimalBinning(max_n_bins=3, name='Azucar')
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
x_transform_bins = optb.transform(imput_transf_cor['Azucar'].values, metric="bins")
x_transform_bins
## array(['(-inf, 44.35)', '(-inf, 44.35)', '(-inf, 44.35)', ...,
##        '(-inf, 44.35)', '(-inf, 44.35)', '(-inf, 44.35)'], dtype=object)
sns.violinplot(x=x_transform_bins,y=varObjCont,palette='viridis')

Podemos probar a realizar el binning de alguna (o todas) las variables continuas que no parezcan relevantes en el anterior análisis en pos de la mejor capacidad predictiva frente al modelo. Recomendable evaluar la salida sin restricciones del árbol y tal vez refinar con un máximo de categorías de salida en su caso. EL violin, revela las ligerísimas diferencias en distribución que se consiguen en este caso tramificando la variable Azucar.

En cualquier caso, es una opción interesante a considerar para la modelización predictiva.

Cuidado con pasarse de categorías!! Como ya sabemos las variables categóricas “consumen” grados de libertad y añaden k-1 parámetros al modelo, aumentando la complejidad del mismo y puediendo incurrir en problemas de robustez o generalización a nuevos datos. Por ello, hay que controlar el crecimiento de categorías de los factores.

Vamos a valorar la capacidad a priori de esta variable.

# ANOVA Azucar tramificada
vinosCompra['Azucar_rec'] = x_transform_bins
model = ols('Beneficio ~ C(Azucar_rec)', data=vinosCompra).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
anova_table
##                      sum_sq      df         F    PR(>F)
## C(Azucar_rec)  2.477983e+05     2.0  2.488539  0.083134
## Residual       2.486906e+08  4995.0       NaN       NaN

No parece que haya una asociación relevante a nivel estadístico.

Son muchas las posibilidades sobre transformaciones de variables de cara a la modelización predictiva. Una estrategia interesante es aplicar múltiples vías de modelización con distintas posibilidades y comparar los resultados finales, decidiendo en base a la capacidad preditiva de la solución y su interpretabilidad.

Por el momento, nos centraremos en la creación de modelos de regresión para estimar el Beneficio de los vinos con Compra = 1 debido a esa mejor en la distribcuión de la respuesta.

Modelos de regresión lineal para la predicción del beneficio del vino

En esta sección se ajustan distintos modelos de regresión lineal para predecir el beneficio de los vinos con compra = 1. En primer lugar, tomamos la partición training (donde ajustamos el modelo) y test (donde probamos su capacidad).

Partición training-test

La filosofía de muchos de los modelos predictivos en python es la de tener por separado la matriz de predictores y el vector de variable objetivo ya que bajo este esquema están programados los procedimientos. Sin embargo, esto no es algo directo cuando existen variables categóricas en el dataset. En general, estas variables no son aceptadas tal cual en los modelos y hay que realizar un paso previo de creación de variables dummy generando lo que se conoce como la matriz de diseño del modelo.

Recordamos aquí lo que ya hemos comentado. Si tengo 4 categorías en un factor y pretendo generar dummies para todas, qué pasa?? Pues que tenemos una combinación lineal exacta de las variables en la matriz y por tanto no se puede invertir (algo necesario para la estimación por mínimos cuadrados de los parámetros del modelo) y todo se va al garete!! Por ello habría que tener en cuenta esto para generar solamente 3 dummies quedando “implícito” el efecto de la cuarta cuando todas ellas tomen el valor 0. Ok…y En este momento nos preguntamos, y donde va el efecto de esta variable? es decir, como se contabiliza en el modelo? Pues precisamente en nuestra querida constante \(\beta_0\) u ordenada en el origne de la ecuación de regresión. Vale, y donde está esta constante…pues depende..si generamos todo “a mano” es decir la extensión en k-1 dummies, tendríamos que añadir, así mismo esta constante (en realidad un vector de todo 1’s en la matriz de diseño).

# Función necesaria
from sklearn.model_selection import train_test_split

# Creamos 4 objetos: predictores para tr y tst y variable objetivo para tr y tst. 
X_train, X_test, y_train, y_test = train_test_split(imputCompra, varObjCont, test_size=0.2, random_state=42)

# Comprobamos dimensiones
print('Training dataset shape:', X_train.shape, y_train.shape)
## Training dataset shape: (3998, 16) (3998,)
print('Testing dataset shape:', X_test.shape, y_test.shape)
## Testing dataset shape: (1000, 16) (1000,)

Vamos a explorar algunas posibilidades para automatizar estas cosas y evitar de momento la necesidad de generar explícitamente y mucho menos a mano, las matrices de diseño.

Algo muy guay que tiene R es la interfaz fórmula, mediante la cual podemos undicar el modelos que queramos ajustar de una forma muy visual. Por ejemplo, como hemos hecho en ANOVA, Beneficio en función de Región o Etiqueta o los efectos que queramos incluir…no tenemos más que poner ‘Beneficio ~ Efecto1 + Efecto2 +..+EfectoN’ y por debajo se estarán generando estas matrices de diseño con sus dummies en caso de necesidad y su constante y de todo.

Tenemos esta posibilidad en python para los modelos de regresión. Generemos primero el data_train con la variable objetivo dentro y posterioremente hacemos uso de la api para fórmulas de statmodels en su variante para regresión por mínimos cuadrados ordinarios.

# Genero el training con la objetivo dentro 
data_train = X_train.join(y_train)
# Importamos la api para fórmulas (en concreto ols para regresión)
from statsmodels.formula.api import ols 

# Ajusto regresión de ejemplo
results = ols('Beneficio ~ Etiqueta + Acidez',data=data_train).fit()
results.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.370
Model: OLS Adj. R-squared: 0.370
Method: Least Squares F-statistic: 469.7
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:09:52 Log-Likelihood: -26401.
No. Observations: 3998 AIC: 5.281e+04
Df Residuals: 3992 BIC: 5.285e+04
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 288.7290 13.962 20.680 0.000 261.356 316.102
Etiqueta[T.M] 138.7840 14.979 9.265 0.000 109.416 168.152
Etiqueta[T.R] 291.5636 14.522 20.077 0.000 263.092 320.036
Etiqueta[T.B] 452.1047 15.051 30.037 0.000 422.596 481.614
Etiqueta[T.MB] 585.2175 20.101 29.115 0.000 545.809 624.626
Acidez -7.0863 3.685 -1.923 0.055 -14.311 0.139
Omnibus: 131.498 Durbin-Watson: 1.991
Prob(Omnibus): 0.000 Jarque-Bera (JB): 143.931
Skew: 0.461 Prob(JB): 5.57e-32
Kurtosis: 3.122 Cond. No. 13.8


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Se peude comprobar que el funcionamiento es adecuado en tanto en cuanto la variable categórica Etiqueta aparece con 5-1 = 4 efectos en el modelo (siendo el restante ‘MM’ recogido por la constante \(\beta_0\)). Además, vemos que se conserva el órden de las categorías que hemos especificado y también que el efetco tiene pinta de monótono creciente con etiqueta puesto que los betas van aumentando con el aumento de etiqueta. Si nos ponemos a valorar, pues es un modelo bastante malo con un R2 de 0,37 pero todos los parámetros son significativos a nivel estadístico (es decir, su parámetro es distinto de 0 o el IC no contiene al valor 0). Ajuste muy pobre con significación paramétrica.

Lo que siempre haremos es un modelo completo de referencia para valorar la capacidad (sobre training) de este y tomarla como base para la generación de modelos más sencillos que mantengan valor predictivo en training y lo mejoren en test.

Modelo completo de referencia

Para utilizar la interfaz fórmula proporcionada por la api, tenemos que especificar todos los efectos…cosa que no me gusta nada (R tiene la posibilidad de poner ‘Beneficio ~ .’ y se sobreentiende que quieres todos los efectos) pero la vida es durilla a veces jeje.

Venga pues vamos a facilitarnos la vida un poco y generamos una función que concatene todos los efectos como string e incluso podamos eliminar algunos a placer, obteniendo la fórmula de manera automática.

# Función para generar la fórmula por larga que sea
def ols_formula(df, dependent_var, *excluded_cols):
    df_columns = list(df.columns.values)
    df_columns.remove(dependent_var)
    for col in excluded_cols:
        df_columns.remove(col)
    return dependent_var + ' ~ ' + ' + '.join(df_columns)

# Aplicamos a fórmula de modelo completo
form=ols_formula(data_train,'Beneficio')
form
## 'Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2'
# Ajusto regresión según fórmula completa
modeloC = ols(form,data=data_train).fit()
modeloC.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.435
Model: OLS Adj. R-squared: 0.431
Method: Least Squares F-statistic: 122.2
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:10:01 Log-Likelihood: -26186.
No. Observations: 3998 AIC: 5.242e+04
Df Residuals: 3972 BIC: 5.259e+04
Df Model: 25
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 273.6801 107.982 2.534 0.011 61.974 485.386
Etiqueta[T.M] 123.1336 14.323 8.597 0.000 95.052 151.215
Etiqueta[T.R] 252.3407 14.069 17.937 0.000 224.758 279.923
Etiqueta[T.B] 383.3082 14.824 25.857 0.000 354.244 412.372
Etiqueta[T.MB] 495.2993 19.721 25.115 0.000 456.635 533.963
Clasificacion[T.*] 17.9532 9.192 1.953 0.051 -0.068 35.974
Clasificacion[T.**] 63.0280 8.857 7.117 0.000 45.664 80.392
Clasificacion[T.***] 127.4193 9.827 12.966 0.000 108.152 146.686
Clasificacion[T.****] 205.0624 13.503 15.187 0.000 178.590 231.535
Region[T.2.0] 12.1224 6.623 1.830 0.067 -0.863 25.107
Region[T.3.0] 0.9415 6.565 0.143 0.886 -11.929 13.812
prop_missings[T.7.6923076923076925] 3.0431 6.504 0.468 0.640 -9.708 15.794
prop_missings[T.16.666666666666664] 21.9144 17.222 1.272 0.203 -11.850 55.679
prop_missings[T.27.27272727272727] -49.3978 64.388 -0.767 0.443 -175.634 76.838
Acidez -7.2167 3.510 -2.056 0.040 -14.099 -0.334
AcidoCitrico -1.6689 3.222 -0.518 0.605 -7.986 4.648
Azucar -0.1400 0.079 -1.778 0.076 -0.294 0.014
CloruroSodico -5.9563 8.484 -0.702 0.483 -22.590 10.678
Densidad -61.1204 106.290 -0.575 0.565 -269.509 147.269
pH 3.8509 4.030 0.956 0.339 -4.050 11.752
Sulfatos 0.7644 2.885 0.265 0.791 -4.892 6.421
Alcohol 4.3407 0.751 5.779 0.000 2.868 5.813
CalifProductor -6.3493 2.502 -2.538 0.011 -11.254 -1.445
PrecioBotella -0.3459 1.844 -0.188 0.851 -3.961 3.270
aleatorio -3.4986 9.380 -0.373 0.709 -21.888 14.891
aleatorio2 12.1311 9.311 1.303 0.193 -6.123 30.386
Omnibus: 113.210 Durbin-Watson: 1.965
Prob(Omnibus): 0.000 Jarque-Bera (JB): 122.531
Skew: 0.429 Prob(JB): 2.47e-27
Kurtosis: 3.029 Cond. No. 1.96e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.96e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Modelo completo bastante sobreparametrizado por lo que se ve, con muchos efectos no significativos cuya aportación al R2 será seguramente despreciable y puedan ser eliminados y con una capacidad bastante baja con R2 de 0.43…mal augurio. Conjunto de datos dificil de modelizar.. Pero hay que tirar palante con lo que tenemos y luchar!

Importacia de las variables

Vamos a tratar de extraer la importancia de las variables según su aportación al R2 total del modelo. Con la función relativeImp podemos realizar esta labor. Sin embargo, de nuevo u pequeño obstáculo…solamente funciona con variables continuas..vaya!

Probamos.

from relativeImp import relativeImp

# Nombres de las variables continuas
names=X_train.select_dtypes(include=np.number).columns.tolist()

# Calcular importancia relativa de las variable continuas en modelo (solo con variables continuas!)
df_results = relativeImp(data_train, outcomeName = 'Beneficio', driverNames = names)

df_results.sort_values(by='normRelaImpt', ascending=False)
##             driver  rawRelaImpt  normRelaImpt
## 7          Alcohol     0.006235     52.287319
## 8   CalifProductor     0.002157     18.086350
## 0           Acidez     0.001874     15.715050
## 11      aleatorio2     0.000492      4.129008
## 2           Azucar     0.000465      3.898723
## 5               pH     0.000180      1.507559
## 9    PrecioBotella     0.000158      1.328349
## 4         Densidad     0.000151      1.266428
## 6         Sulfatos     0.000074      0.623407
## 3    CloruroSodico     0.000065      0.543347
## 10       aleatorio     0.000037      0.313968
## 1     AcidoCitrico     0.000036      0.300494

Con esta funcionalidad siempre podemos calcular la importacia relativa según la aportación al R2 global de cada una de las variables de un modelo especifico. Sin embargo, la función no acepta variables categóricas por lo que, para poder tener información sobre la importacia del modelo completo, es necesario realizar la extensión en dummies (algo que sucede en el interior del modelo de regresión con la api de fórmula).

En el modelo solamente con variables continuas (con un R2 lamentable como de 0.01) las más importantes son Alcohol, CalifProductor y Acidez.

Lo que nos gustaría es poder ver esto para el modelo completo incluyendo las nominales. Para ello, no nos queda otra alternativa que generar de manera explícita las matrices de diseño del modelo.

El paquete pasty resulta muy útil para generar las matrices de diseño de un modelo basadas en una fórmula concreta. Esto nos permitirá utilizar cualquier modelo de predicción de python que no acepte directamente la interfaz fórmula. Es evidente que podríamos hacer esta operación “a mano”, generando la extensión en dummies de los factores y la posterior adición de una constante \(\beta_0\) para que se refleje en el modelo, pero pienso que es más intuitivo utilizar la fórmula y así evitar errores.

import patsy

# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices(form, data_train, return_type='dataframe')

# Ahora podemos aplicar la función "oficial" de statmodels OLS (con formato y,X)
model=sm.OLS(y,X).fit()
model.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.435
Model: OLS Adj. R-squared: 0.431
Method: Least Squares F-statistic: 122.2
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:10:09 Log-Likelihood: -26186.
No. Observations: 3998 AIC: 5.242e+04
Df Residuals: 3972 BIC: 5.259e+04
Df Model: 25
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 273.6801 107.982 2.534 0.011 61.974 485.386
Etiqueta[T.M] 123.1336 14.323 8.597 0.000 95.052 151.215
Etiqueta[T.R] 252.3407 14.069 17.937 0.000 224.758 279.923
Etiqueta[T.B] 383.3082 14.824 25.857 0.000 354.244 412.372
Etiqueta[T.MB] 495.2993 19.721 25.115 0.000 456.635 533.963
Clasificacion[T.*] 17.9532 9.192 1.953 0.051 -0.068 35.974
Clasificacion[T.**] 63.0280 8.857 7.117 0.000 45.664 80.392
Clasificacion[T.***] 127.4193 9.827 12.966 0.000 108.152 146.686
Clasificacion[T.****] 205.0624 13.503 15.187 0.000 178.590 231.535
Region[T.2.0] 12.1224 6.623 1.830 0.067 -0.863 25.107
Region[T.3.0] 0.9415 6.565 0.143 0.886 -11.929 13.812
prop_missings[T.7.6923076923076925] 3.0431 6.504 0.468 0.640 -9.708 15.794
prop_missings[T.16.666666666666664] 21.9144 17.222 1.272 0.203 -11.850 55.679
prop_missings[T.27.27272727272727] -49.3978 64.388 -0.767 0.443 -175.634 76.838
Acidez -7.2167 3.510 -2.056 0.040 -14.099 -0.334
AcidoCitrico -1.6689 3.222 -0.518 0.605 -7.986 4.648
Azucar -0.1400 0.079 -1.778 0.076 -0.294 0.014
CloruroSodico -5.9563 8.484 -0.702 0.483 -22.590 10.678
Densidad -61.1204 106.290 -0.575 0.565 -269.509 147.269
pH 3.8509 4.030 0.956 0.339 -4.050 11.752
Sulfatos 0.7644 2.885 0.265 0.791 -4.892 6.421
Alcohol 4.3407 0.751 5.779 0.000 2.868 5.813
CalifProductor -6.3493 2.502 -2.538 0.011 -11.254 -1.445
PrecioBotella -0.3459 1.844 -0.188 0.851 -3.961 3.270
aleatorio -3.4986 9.380 -0.373 0.709 -21.888 14.891
aleatorio2 12.1311 9.311 1.303 0.193 -6.123 30.386
Omnibus: 113.210 Durbin-Watson: 1.965
Prob(Omnibus): 0.000 Jarque-Bera (JB): 122.531
Skew: 0.429 Prob(JB): 2.47e-27
Kurtosis: 3.029 Cond. No. 1.96e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.96e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
# Ajusto regresión de ejemplo
results = ols('Beneficio ~ aleatorio2 + Etiqueta',data=data_train).fit()
results.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.370
Model: OLS Adj. R-squared: 0.369
Method: Least Squares F-statistic: 469.3
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:10:14 Log-Likelihood: -26402.
No. Observations: 3998 AIC: 5.282e+04
Df Residuals: 3992 BIC: 5.285e+04
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 278.1061 14.725 18.887 0.000 249.237 306.976
Etiqueta[T.M] 139.7310 14.976 9.330 0.000 110.369 169.093
Etiqueta[T.R] 292.7666 14.514 20.172 0.000 264.312 321.222
Etiqueta[T.B] 453.0996 15.044 30.118 0.000 423.605 482.595
Etiqueta[T.MB] 586.8203 20.095 29.202 0.000 547.422 626.218
aleatorio2 14.7823 9.784 1.511 0.131 -4.401 33.965
Omnibus: 131.244 Durbin-Watson: 1.993
Prob(Omnibus): 0.000 Jarque-Bera (JB): 143.626
Skew: 0.460 Prob(JB): 6.49e-32
Kurtosis: 3.122 Cond. No. 14.3


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Podemos comprobar que el modelo es exactamente el mismo que el ajustado mediante la fórmula y la función ols, pero ahora resulta que las matrices de diseño se especifican explicitamente con lo que podemos utilizar la importancia relativa de variables ya que no existen factores. Veamos..

# Nombres de predictores (en modo dummy) donde quitamos la constante
names=X.columns.tolist()[1:]

# Calculamos importancia relativa
df_results = relativeImp(X.join(y), outcomeName = 'Beneficio', driverNames = names)

# Ordenamos valores 
df_results.sort_values(by='normRelaImpt', ascending=False)
##                                  driver  rawRelaImpt  normRelaImpt
## 2                         Etiqueta[T.B]     0.141921     32.645204
## 3                        Etiqueta[T.MB]     0.072831     16.752895
## 7                 Clasificacion[T.****]     0.057918     13.322638
## 0                         Etiqueta[T.M]     0.048181     11.082912
## 6                  Clasificacion[T.***]     0.046103     10.604822
## 1                         Etiqueta[T.R]     0.034920      8.032400
## 4                    Clasificacion[T.*]     0.014820      3.408889
## 5                   Clasificacion[T.**]     0.006730      1.548010
## 20                              Alcohol     0.005406      1.243454
## 21                       CalifProductor     0.001583      0.364180
## 13                               Acidez     0.001324      0.304557
## 8                         Region[T.2.0]     0.000728      0.167488
## 15                               Azucar     0.000522      0.119964
## 11  prop_missings[T.16.666666666666664]     0.000413      0.095110
## 24                           aleatorio2     0.000349      0.080353
## 12   prop_missings[T.27.27272727272727]     0.000261      0.059923
## 18                                   pH     0.000213      0.048928
## 9                         Region[T.3.0]     0.000111      0.025580
## 16                        CloruroSodico     0.000104      0.023823
## 17                             Densidad     0.000099      0.022765
## 14                         AcidoCitrico     0.000056      0.012985
## 22                        PrecioBotella     0.000046      0.010545
## 19                             Sulfatos     0.000045      0.010264
## 23                            aleatorio     0.000034      0.007854
## 10  prop_missings[T.7.6923076923076925]     0.000019      0.004455
# Gráfico de importancia relativa en base al R2
px.bar(df_results,x='normRelaImpt',y='driver',title='Importancia relativa por aportación al R2').update_yaxes(categoryorder="total ascending").show()

Ahora si que si! Obtenemos un gráfico de importancia de las variables del modelo completo con el efecto de las categorías de las nominales por separado. Claramente lo que se podía esperar, las categorías de Clasificación y Etiqueta son las de mayor importancia para el modelo.

Vamos a generar un proceso backward eliminando variables secuencialmente según el p-valor.

# Proceso backward de eliminación de efectos según p-valor
form2=ols_formula(data_train,'Beneficio','prop_missings','aleatorio','aleatorio2',
                  'PrecioBotella','Sulfatos','AcidoCitrico','Densidad',
                  'CloruroSodico','pH','Region','Azucar')

# Ajusto regresión sin prop_missings
modeloC2 = ols(form2,data=data_train).fit()
modeloC2.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.433
Model: OLS Adj. R-squared: 0.431
Method: Least Squares F-statistic: 276.5
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:10:25 Log-Likelihood: -26193.
No. Observations: 3998 AIC: 5.241e+04
Df Residuals: 3986 BIC: 5.248e+04
Df Model: 11
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 232.5770 17.809 13.059 0.000 197.661 267.493
Etiqueta[T.M] 123.9191 14.299 8.666 0.000 95.885 151.954
Etiqueta[T.R] 253.0856 14.042 18.023 0.000 225.555 280.616
Etiqueta[T.B] 384.7283 14.798 26.000 0.000 355.717 413.740
Etiqueta[T.MB] 496.4143 19.690 25.211 0.000 457.810 535.019
Clasificacion[T.*] 17.7085 9.173 1.930 0.054 -0.276 35.693
Clasificacion[T.**] 62.7754 8.842 7.100 0.000 45.440 80.111
Clasificacion[T.***] 126.8453 9.814 12.925 0.000 107.604 146.087
Clasificacion[T.****] 204.3448 13.491 15.147 0.000 177.895 230.794
Acidez -6.9836 3.501 -1.995 0.046 -13.848 -0.119
Alcohol 4.4198 0.750 5.895 0.000 2.950 5.890
CalifProductor -6.8412 2.481 -2.757 0.006 -11.705 -1.977
Omnibus: 114.292 Durbin-Watson: 1.967
Prob(Omnibus): 0.000 Jarque-Bera (JB): 123.837
Skew: 0.431 Prob(JB): 1.29e-27
Kurtosis: 3.016 Cond. No. 135.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
from sklearn.metrics import mean_squared_error, r2_score

# Predicciones para test modelo completo
vinos_y_pred = modeloC.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, vinos_y_pred)))
# The coefficient of determination: 1 is perfect prediction
## Mean squared error: 161.55
print("Coefficient of determination: %.2f" % r2_score(y_test, vinos_y_pred))

# Gráfico de real frente a prediccion para alguna variable
#plt.scatter(X_test.Acidez, y_test, color="black")
#plt.scatter(X_test.Acidez, vinos_y_pred, color="blue")
## Coefficient of determination: 0.44
from sklearn.metrics import mean_squared_error, r2_score

# Predicciones para test modelo C2
vinos_y_pred = modeloC2.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, vinos_y_pred)))
# The coefficient of determination: 1 is perfect prediction
## Mean squared error: 161.31
print("Coefficient of determination: %.2f" % r2_score(y_test, vinos_y_pred))

# Gráfico de real frente a prediccion para alguna variable
#plt.scatter(X_test.Acidez, y_test, color="black")
#plt.scatter(X_test.Acidez, vinos_y_pred, color="blue")
## Coefficient of determination: 0.44

Ajuste por validación cruzada repetida

En esta sección vamos a valorar el ajuste de modelos por el proceso de validación cruzada repetida, en el que se generarán n_splits particiones del archivo, repitiendo el proceso con distintas semillas n_repeats veces. De esta forma, obtenemos n_splits x n_repeats modelos por fórmula especificada, pudiendo así promediar los resultados obtenidos bajo muchas muestras de training-test.

En este esquema, todas las observaciones del archivo son utilizadas en algún training para el ajuste del modelo yen algún test para su evaluación pero nunca simultáneamente. Así el modelo tiene instancias cambiantes para el ajuste de parámetros y para la predicción.

# from sklearn.model_selection import cross_val_score
# from sklearn.model_selection import RepeatedKFold
# from sklearn.linear_model import LinearRegression

model = LinearRegression()
#model = sm.OLS(y,X)

# Establecemos esquema de validación fijando random_state (reproducibilidad)
cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=12345)

# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X, y, cv=cv)

# Sesgo y varianza
print('Coeficioente de detrminación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
## Coeficioente de detrminación R2: 0.425 (0.023)
sns.violinplot(y=scores,palette='viridis')

Modelo final backward

# Generamos las matrices de diseño según la fórmula de modelo completo
y2, X2 = patsy.dmatrices(form2, vinosCompra, return_type='dataframe')

# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X2, y2, cv=cv)

# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
## Coeficioente de detrminación R2 Modelo Final Backward: 0.430 (0.020)
sns.violinplot(y=scores,palette='viridis')

Modelo con interacción Etiqueta-Clasificación

Siendo variables que miden de alguna forma la calidad del vino, es posible que exista cierta interacción entre ellas, es decir que la pertenencia conjunta a determinada etiqueta y clasificación tenga efectos diferenciadores en el beneficio que produce el vino.

form3 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion'

# Ajusto regresión sin prop_missings
modeloC3 = ols(form3,data=data_train).fit()
modeloC3.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.439
Model: OLS Adj. R-squared: 0.435
Method: Least Squares F-statistic: 119.6
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:10:51 Log-Likelihood: -26170.
No. Observations: 3998 AIC: 5.239e+04
Df Residuals: 3971 BIC: 5.256e+04
Df Model: 26
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 214.0458 25.322 8.453 0.000 164.401 263.690
Etiqueta[T.M] 123.3191 25.968 4.749 0.000 72.408 174.231
Etiqueta[T.R] 265.0898 26.056 10.174 0.000 214.006 316.174
Etiqueta[T.B] 481.7667 33.008 14.596 0.000 417.053 546.481
Etiqueta[T.MB] 694.7219 60.959 11.396 0.000 575.207 814.237
Clasificacion[T.*] 48.5567 30.214 1.607 0.108 -10.680 107.793
Clasificacion[T.**] 91.3712 40.012 2.284 0.022 12.926 169.816
Clasificacion[T.***] 95.6293 60.987 1.568 0.117 -23.940 215.198
Clasificacion[T.****] 124.0220 17.909 6.925 0.000 88.911 159.133
Etiqueta[T.M]:Clasificacion[T.*] -8.2539 33.620 -0.246 0.806 -74.168 57.660
Etiqueta[T.R]:Clasificacion[T.*] -22.9217 33.481 -0.685 0.494 -88.563 42.720
Etiqueta[T.B]:Clasificacion[T.*] -144.7659 40.925 -3.537 0.000 -225.003 -64.529
Etiqueta[T.MB]:Clasificacion[T.*] -112.5326 87.554 -1.285 0.199 -284.189 59.123
Etiqueta[T.M]:Clasificacion[T.**] -13.2018 42.632 -0.310 0.757 -96.785 70.381
Etiqueta[T.R]:Clasificacion[T.**] -19.0431 42.242 -0.451 0.652 -101.861 63.774
Etiqueta[T.B]:Clasificacion[T.**] -113.6291 47.330 -2.401 0.016 -206.423 -20.835
Etiqueta[T.MB]:Clasificacion[T.**] -224.5732 74.642 -3.009 0.003 -370.914 -78.233
Etiqueta[T.M]:Clasificacion[T.***] 74.5139 64.540 1.155 0.248 -52.021 201.049
Etiqueta[T.R]:Clasificacion[T.***] 26.5701 62.673 0.424 0.672 -96.304 149.444
Etiqueta[T.B]:Clasificacion[T.***] -38.5356 66.105 -0.583 0.560 -168.139 91.068
Etiqueta[T.MB]:Clasificacion[T.***] -161.0492 85.715 -1.879 0.060 -329.099 7.000
Etiqueta[T.M]:Clasificacion[T.****] 160.8211 44.287 3.631 0.000 73.994 247.648
Etiqueta[T.R]:Clasificacion[T.****] 101.9520 25.367 4.019 0.000 52.218 151.686
Etiqueta[T.B]:Clasificacion[T.****] -4.1577 27.740 -0.150 0.881 -58.544 50.229
Etiqueta[T.MB]:Clasificacion[T.****] -134.5934 52.123 -2.582 0.010 -236.784 -32.403
Acidez -7.2240 3.500 -2.064 0.039 -14.086 -0.362
Alcohol 4.4645 0.748 5.967 0.000 2.998 5.931
CalifProductor -6.5027 2.477 -2.626 0.009 -11.358 -1.647
Omnibus: 108.268 Durbin-Watson: 1.966
Prob(Omnibus): 0.000 Jarque-Bera (JB): 116.838
Skew: 0.419 Prob(JB): 4.26e-26
Kurtosis: 2.992 Cond. No. 1.05e+16


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.82e-27. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
# Generamos las matrices de diseño según la fórmula de modelo completo
y3, X3 = patsy.dmatrices(form3, vinosCompra, return_type='dataframe')

# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X3, y3, cv=cv)

# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
## Coeficioente de detrminación R2 Modelo Final Backward: 0.422 (0.027)
sns.violinplot(y=scores,palette='viridis')

Llegados a este punto se observa que, si bien el modelo con la interacción es bastante bueno, tiene un cruce de categorías de la interacción que no contiene observaciones en data_train para poder ser estimado su parámetro…esto es un fallo en la especificación del modelo y será difícil mantenerlo tal cual puesto que a nivel interpretativo carece de la robustez necesaria.

Una estrategia que se puede seguir es unir alguna de las categorías implicadas en esta interacción peligrosa. En este caso podemos unir las categorías de clasificación más elevadas (se puede probar trabajando sobre las etiquetas!!).

Cuando hacemos lago así, hay que tener la precaución de hacer el cambio en el dataset completo y volver a generar la partición para que los datos se actualicen en training y también en test.

Actualizaremos el último modelo con esta unión de categorías.

imputCompra['Etiqueta_r'] = vinosCompra.Etiqueta_r
imputCompra['Azucar_rec'] = vinosCompra.Azucar_rec

# Creamos 4 objetos: predictores para tr y tst y variable objetivo para tr y tst. 
X_train, X_test, y_train, y_test = train_test_split(imputCompra, varObjCont, test_size=0.2, random_state=42)

# Genero el training con la objetivo dentro 
data_train_r = X_train.join(y_train)

form4 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion'

# Generamos las matrices de diseño según la fórmula de modelo completo
y4, X4 = patsy.dmatrices(form4, vinosCompra, return_type='dataframe')

# Obtenemos los resultados de R2 para cada partición tr-tst
scores = cross_val_score(model, X4, y4, cv=cv)

# Sesgo y varianza
print('Coeficioente de detrminación R2 Modelo Final Backward interacción Rec: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
## Coeficioente de detrminación R2 Modelo Final Backward interacción Rec: 0.421 (0.019)
sns.violinplot(y=scores,palette='viridis')

Comparación por validación cruzada

Una vez explorados varios modelos manuales, es bueno comprobar sus capacidades en un esquema de validación cruzada repetida con la intención de comprobar su estabilidad ante el remuestreo.

El objetivo fundamental en esta parte es seleccionar el mejor modelo en relación sesgo-varianza y complejidad.

Vamos a generar una función que automatice un poco el proceso de comparación por validación cruzada, de tal forma que pueda aplicar sobre cualquier fórmula epecificada. Luego será cuestión de enlistar todas las fómulas que queramos y pasar la función con un map().

# Función para comparación por validación cruzada
def cross_val(formula, data, seed=12345):
      # Generamos las matrices de diseño según la fórmula de modelo completo
      y, X = patsy.dmatrices(formula, data, return_type='dataframe')
  
      # Establecemos esquema de validación fijando random_state (reproducibilidad)
      cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=seed)
  
      # Obtenemos los resultados de R2 para cada partición tr-tst
      scores = cross_val_score(model, X, y, cv=cv)
  
      # Sesgo y varianza
      print('Modelo: ' + formula)
      print('Coeficiente de determinación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
      
      #sns.violinplot(y=scores,palette='viridis')
      
      return(scores)
form5 = 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec'
form6 = 'Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion'
form7 = 'Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec'
form8 = 'Beneficio ~ Etiqueta + Clasificacion + Azucar_rec'


# Creamos lista de fórmulas   
list_form = [form,form2,form3,form4,form5,form6,form7,form8]
list_form

# Aplicamos a toda la lista la función creada (devuelve un dataframe pero está transpuesto)
## ['Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion', 'Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec', 'Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion', 'Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec', 'Beneficio ~ Etiqueta + Clasificacion + Azucar_rec']
list_res = pd.DataFrame(map(lambda x: cross_val(x,vinosCompra, seed=2022),list_form))

# Trasnponer dataframe y pasar de wide a long (creando un factor variable con el nombre de cada fórmula de la lista[0,1,2,3])
## Modelo: Beneficio ~ Acidez + AcidoCitrico + Azucar + CloruroSodico + Densidad + pH + Sulfatos + Alcohol + CalifProductor + PrecioBotella + Etiqueta + Clasificacion + Region + prop_missings + aleatorio + aleatorio2
## Coeficiente de determinación R2: 0.429 (0.019)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion
## Coeficiente de determinación R2: 0.431 (0.018)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta*Clasificacion
## Coeficiente de determinación R2: 0.426 (0.018)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta_r*Clasificacion
## Coeficiente de determinación R2: 0.421 (0.019)
## Modelo: Beneficio ~ Acidez + Alcohol + CalifProductor + Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.431 (0.018)
## Modelo: Beneficio ~ Alcohol + CalifProductor + Etiqueta + Clasificacion
## Coeficiente de determinación R2: 0.430 (0.018)
## Modelo: Beneficio ~ CalifProductor + Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.426 (0.018)
## Modelo: Beneficio ~ Etiqueta + Clasificacion + Azucar_rec
## Coeficiente de determinación R2: 0.425 (0.019)
results = list_res.T.melt()
results.columns = ['Modelo','R2']
results.head()
##    Modelo        R2
## 0       0  0.394846
## 1       0  0.453454
## 2       0  0.388224
## 3       0  0.427713
## 4       0  0.457318

Una vez generados los 100 modelos por fórmula con distintas particiones de trainig test bajo un esquema de validación cruzada repetida 5 Folds 20 Repeats, juntamos los resultados de los “resamples” en un dataset para generar el gráfico de boxplots de las distribuciones de R2 a través de los modelos y valorar la relación sesgo-varianza.

# Boxplot paralelo para comparar
sns.boxplot(x='Modelo',y='R2',data=results,palette='viridis')

En este punto hay que decidir un modelo final.

Todo indica que el mejor modelo en todos los sentidos es el modelo final backwar sin interacción. Modelo sencillo y con capacidad predictiva mayor en términos generales que sus competidores que además son más complejos. Claramente nos quedamos con el indicado como 1 que es el modelo de fórmula form2.

Ajuste del modelo final. Interpretación de parámetros

Una vez escogido el modelo final. Hay que interpretar sus parámetros teniendo en cuenta que la estimación más robusta de los mismos será la que se obtiene utilizando el datset completo y no solamente las instancias del conjunto de training.

De esta forma, se ajusta el modelo final a los datos completos.

# Ajusto regresión sin prop_missings
modelo_final = ols(form2,data=vinosCompra).fit()
modelo_final.summary()
OLS Regression Results
Dep. Variable: Beneficio R-squared: 0.435
Model: OLS Adj. R-squared: 0.433
Method: Least Squares F-statistic: 348.5
Date: do., 30 oct. 2022 Prob (F-statistic): 0.00
Time: 21:11:23 Log-Likelihood: -32696.
No. Observations: 4998 AIC: 6.542e+04
Df Residuals: 4986 BIC: 6.549e+04
Df Model: 11
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 237.4493 15.767 15.060 0.000 206.539 268.359
Etiqueta[T.M] 123.4425 12.656 9.753 0.000 98.631 148.254
Etiqueta[T.R] 255.8916 12.431 20.585 0.000 231.522 280.261
Etiqueta[T.B] 385.2418 13.117 29.370 0.000 359.527 410.957
Etiqueta[T.MB] 494.7721 17.662 28.014 0.000 460.147 529.397
Clasificacion[T.*] 20.7499 8.105 2.560 0.010 4.860 36.640
Clasificacion[T.**] 66.6822 7.829 8.518 0.000 51.334 82.030
Clasificacion[T.***] 124.9397 8.658 14.431 0.000 107.967 141.913
Clasificacion[T.****] 208.6767 11.971 17.433 0.000 185.209 232.144
Acidez -7.3930 3.107 -2.379 0.017 -13.485 -1.301
Alcohol 4.1374 0.665 6.218 0.000 2.833 5.442
CalifProductor -8.4882 2.216 -3.831 0.000 -12.832 -4.145
Omnibus: 133.624 Durbin-Watson: 2.006
Prob(Omnibus): 0.000 Jarque-Bera (JB): 144.045
Skew: 0.415 Prob(JB): 5.26e-32
Kurtosis: 3.045 Cond. No. 136.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Conclusiones que podemos sacar para aquellos vinos con Compra = 1, es decir beneficio distinto de 0:

  1. Los vinos con Clasificación ** ** tienen un beneficio estimado medio 208.67 unidades superior a los vinos de clasificación Desconocida. El rango de variación esperado es [185.209, 232.144].

  2. El beneficio estimado medio de los vinos de Clasificación * es 20.74 unidades superior al de los vinos con Clasificación Deconocida. El rango de variación esperado es [4.860, 36.640].

  3. El beneficio estimado medio de los vinos de Etiqueta MB, es 494.77 unidades mayor que el de los vinos con etiqueta MM. El rango de variación esperado es [460.147, 529.397].

  4. El aumento unitario de la Calificación del Productor, produce una disminución del beneficio estimado del vino de 8.488 unidades. Puede resutar sorprendente pero hay que tener en cuenta que la variación es a un nivel marginal, con lo cual se valora qué pasa con el aumento unitario de la calificación del propio productor para vinos con la misma clasificación, etiqueta y niveles de alcohol y acidez.

  5. El aumento unitario en el nivel de Acidez produce una *disminución del Beneficio esperado de 7.39 unidades**

  6. El aumento unitario en el nivel de Alcohol produce un *aumento del Beneficio esperado de 4.13 unidades**

Todas estas afirmaciones se hacen a constancia de todas las demás variables involucradas en el modelo (ceteris paribus), es decir, se valora el cambio marginal que tiene cada variable para valores fijos de todas las demás.

Así, ese aumento en los vinos de **** se da en aquellos que tienen la misma etiqueta, misma calificación del productor y mismos niveles de acidez y alcohol.